%% This code runs the Monte Carlo simulations from varying Rscale and cost
% plot a subsample of the simulation results.
close all; clear;
clc;

% define saved files from running this code.
fig1_save = 'fig5a_supply_quantiles.png';
fig2_save = 'fig5b_backP_hist.png';

%% load data and estimates from est1 & est2
load est1_output_2006.mat
load est2_output1_c.mat

% set model parameters
run set_parameters.m

%% distribution of Rscale
Rscale0 = R1_scale; % middle value
RscaleL = 0.5*R1_scale; % low
RscaleH = 1.5*R1_scale; % high
% RscaleL = 0.5; RscaleH = 4; % asymmetric distribution
pd_Rscale = makedist('Triangular',RscaleL,Rscale0,RscaleH);

%% distribution of cost
cost = c; se_cost = se_c;
costL = cost-2*se_cost; % low
costH = cost+2*se_cost; % high
pd_cost = makedist('Triangular',costL,cost,costH);

%% prepare for Monte Carlo simulation
Imax = 500; % number of MC iterations (500)
M = 1000; % number of random draws (1000)

%% plot a sample of simulation (1 MC iteration)
fig1 = figure('Position',[0 0 700 600]);
hold on
% value of random parameter(s) ***********
rng('default') % for replication
Rscale_mc = random(pd_Rscale,M,1);
cost_mc = random(pd_cost,M,1);
backP = zeros(M,2);
Supply = cell(M,1); % collection of supply curves

for m = 1:M
    if mod(m,50) == 0, fprintf('.')
        end
    Rscale = Rscale_mc(m);
    cc = cost_mc(m);
    [H_sim,P_sim,~,backP_mc] = sim_supply(Rscale,q0,Theta,cc);
    plot((1e-3)*H_sim,P_sim,'Color',0.6*[1 1 1]);
    backP(m,:) = [m backP_mc];
    Supply{m} = [H_sim P_sim]; % store simulated supply
end
ylim([0 50]); xlim([0 Inf])
hold off

% quantiles of backP
backP_q0 = quantile(backP(:,2),0.5);
backP_q1 = quantile(backP(:,2),0.025);
backP_q2 = quantile(backP(:,2),0.975);

[~, indx_q0] = min((backP(:,2) - backP_q0).^2);
[~, indx_q1] = min((backP(:,2) - backP_q1).^2);
[~, indx_q2] = min((backP(:,2) - backP_q2).^2);

supply_q0 = Supply{backP(indx_q0,1)};
supply_q1 = Supply{backP(indx_q1,1)};
supply_q2 = Supply{backP(indx_q2,1)};

supply_q0(supply_q0 <=0) = nan;
supply_q1(supply_q1 <=0) = nan;
supply_q2(supply_q2 <=0) = nan;

figure(fig1)
hold on
h0 = plot((1e-3)*supply_q0(:,1),supply_q0(:,2),'k','LineWidth',3);
h1 = plot((1e-3)*supply_q1(:,1),supply_q1(:,2),'k-.','LineWidth',3);
h2 = plot((1e-3)*supply_q2(:,1),supply_q2(:,2),'k-.','LineWidth',3);
lgd = legend([h0 h1],{'median','2.5 & 97.5 quantile'},...
    'Location','best','FontSize',12);
xlabel 'Catch Quantity (metric tons)'
ylabel 'Price ($/kg)'
title '(a) a simulated sample of steady-state supply of fish'
hold off

saveas(fig1,fig1_save)

%% MC simulation (Imax iterations)
backP = zeros(M,Imax); % sample of backward bending prices
rng('default') % for replication
for iter = 1:Imax
    fprintf('\nMC iter (%g)', iter)
    % value of random parameter(s) ***********
    Rscale_mc = random(pd_Rscale,M,1);
    cost_mc = random(pd_cost,M,1);
    
    for m = 1:M
        if mod(m,50) == 0, fprintf('.')
        end
        Rscale = Rscale_mc(m);
        cc = cost_mc(m);
        [H_sim,P_sim,Hmax,backP_mc] = sim_supply(Rscale,q0,Theta,cc);
        backP(m,iter) = backP_mc;
    end
end

% Quantiles of backward bending Price 
backP_q0 = quantile(backP,0.5,1);
backP_q1 = quantile(backP,0.025,1);
backP_q2 = quantile(backP,0.975,1);

%% plot histogram of backP quantiles 
fig2 = figure('Position',[50 50 700 400]);
hold on
h1 = histogram(backP_q0,20,'FaceColor',[0 0.4 0.7],'EdgeColor',[0 0.4 0.7]);
h2 = histogram(backP_q1,20,'FaceColor',[0.4 0.7 0],'EdgeColor',[0.4 0.7 0]);
h3 = histogram(backP_q2,20,'FaceColor',[0.4 0 0.7],'EdgeColor',[0.4 0 0.7]);
hold on
str1 = sprintf('Median\n %.2f',mean(backP_q0)); %gtext(str1,'Color',[0 0.4 0.7])
str2 = sprintf('2.5-Quantile\n %.2f',mean(backP_q1)); %gtext(str2,'Color',[0.4 0.7 0])
str3 = sprintf('97.5-Quantile\n %.2f',mean(backP_q2));%gtext(str3,'Color',[0.4 0 0.7])
title '(b) Histogram of price quantiles and their sample means'
xlabel 'price ($/kg)'
ylim([0 80])
hold off
saveas(fig2,fig2_save)